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The time evolution of a large- amplitude electromagnetic (EM) wave injected vertically into the 
overhead ionosphere is studied numerically. The EM wave has a carrier frequency of 5 MHz and 
is modulated as a Gaussian pulse with a width of approximately 0.1 milliseconds and a vacuum 
amplitude of 1.5 V/m at 50 km. This is a fair representation of a modulated radio wave transmitted 
from a typical high-power HF broadcast station on the ground. The pulse is propagated through 
the neutral atmosphere to the critical points of the ionosphere, where the L-O and R-X modes 
are reflected, and back to the neutral atmosphere. We observe mode conversion of the L-O mode 
to electrostatic waves, as well as harmonic generation at the turning points of both the R-X and 
L-O modes, where their amplitudes rise to several times the original ones. The study has relevance 
for ionospheric interaction experiments in combination with ground-based and satellite or rocket 
observations. 



I. INTRODUCTION 



Pulsed high-frequency (HF) electromagnetic (EM) 
waves from transmitters on the ground are regularly used 
for sounding the density profile and drift velocity of the 
overehead ionosphere [Hunsucker, 1991; Reinisch et ai., 
1995, Reinisch, 1996]. In 1971, it was shown theoreti- 
cally by Perkins and Raw [1971] that if the injected HF 
radio beams are strong enough, weak-turbulence para- 
metric instabilities in the ionospheric plasma of the type 
predicted by Silin [1965] and DuBois and Goldman [1965] 
would be excited. Ionospheric modification experiments 
by a high-power HF radio wave at Platteville in Colorado 
[Utlaut, 1970], using ionosonde recordings and photomet- 
ric measurements of artificial airglow, demonstrated the 
heating of electrons, the deformation in the traces on 
ionosonde records, the excitation of spread F, etc., af- 
ter the HF transmitter was turned on. The triggering 
of weak-turbulence parametric instabilities in the iono- 
sphere was first observed in 1970 in experiments on the 
interaction between powerful HF radio beams and the 
ionospheric plasma, conducted at Arecibo, Puerto Rico, 
using a scatter radar diagnostic technique [Wong and 
Taylor, 1971; Carlson et ai, 1972]. A decade later it 
was found experimentally in Troms0 that, under simi- 
lar experimental conditions as in Arecibo, strong, sys- 
tematic, structured, wide-band secondary HF radiation 
escapes from the interaction region [Thide et ai, 1982]. 
This and other observations demonstrated that complex 
interactions, including weak and strong EM turbulence, 
[Leyser, 2001; Thide et at, 2005] and harmonic gener- 
ation [Derblom et ai, 1989; Blagoveshchenskaya et ai, 
1998] are excited in these experiments. 

Numerical simulations have become an important tool 
to understand the complex behavior of plasma turbu- 



lence. Examples include analytical and numerical studies 
of Langmuir turbulence [Robinson, 1997], and of upper- 
hybrid/lower- hybrid turbulence in magnetized plasmas 
[Goodman et ai, 1994; Xi, 2004]. In this Letter, we 
present a full-scale simulation study of the propagation 
of an HF EM wave into the ionosphere, with ionospheric 
parameters typical for the high-latitude EISCAT Heat- 
ing facility in Troms0, Norway. To our knowledge, this 
is the first simulation involving realistic scale sizes of the 
ionosphere and the wavelength of the EM waves. Our 
results suggest that such simulations, which are possible 
with today's computers, will become a powerful tool to 
study HF-induced ionospheric turbulence and secondary 
radiation on a quantitative level for direct comparison 
with experimental data. 



II. MATHEMATICAL MODEL AND 
NUMERICAL SETUP 

We use the MKS system (SI units) in the mathemati- 
cal expressions throughout the manuscript, unless other- 
wise stated. We assume a vertically stratified ion number 
density profile n^z) with a constant geomagnetic field 
Bo directed obliquely to the density gradient. The EM 
wave is injected vertically into the ionosphere, with spa- 
tial variations only in the z direction. Our simple one- 
dimensional model neglects the EM field 1/r falloff (r is 
the distance from the transmitter), the Fresnel pattern 
created obliquely to the z direction by the incident and 
reflected wave, and the the influence on the radio wave 
propagation due to field aligned irregularities in the iono- 
sphere. For the EM wave, the Maxwell equations give 
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where the electron fluid velocity is obtained from the mo- 
mentum equation 

^ = -^--[E + v e x(B + B 1 )] (3) 
at oz m e 

and the electron density is obtained from the Poisson 
equation n c = riio(z) — (eo/e)dE z /dz. Here, z is the 
unit vector in the z direction, c is the speed of light in 
vacuum, e is the magnitude of the electron charge, £o is 
the vacuum permittivity, and m is the electron mass. 
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FIG. 1: The ion density profile, and the electric and magnetic 
field components at time t — ms. 

The number density profile of the immobile ions, 
n i0 (z) = 0.5 x 10 12 exp[-(z - 300) 2 /10 3 ] (z in kilome- 
ters) is shown in the leftmost panel of Fig. [1] Instead of 
modeling a transmitting antenna via a time-dependent 
boundary condition at z — km, we assume that the 
EM pulse has reached the altitude z = 50 km when we 
start our simulation, and we give the pulse as an initial 
condition at time t = s. In the initial condition, we use 
a linearly polarized EM pulse where the carrier wave has 
the wavelength A = 60m (wavenumber k = 0.1047m~ 1 ) 
corresponding to a carrier frequency of /o = 5 MHz (ljq = 
31 x 10 6 s _1 ). The EM pulse is amplitude modulated in 
the form of a Gaussian pulse with a maximum amplitude 
of 1.5 V/m, with the ^-component of the electric field 
set to E x = 1.5exp[-(z - 50) 2 /10 2 ] sin(0.1047 x 10 3 z) 
[z in kilometers) and the y component of the magnetic 
field set to B y = E x /c at t = 0. The other electric and 
magnetic field components are set to zero; see Fig. [T] 
The spatial width of the pulse is approximately 30 km, 
corresponding to a temporal width of 0.1 milliseconds 
as the pulse propagates with the speed of light in the 
neutral atmosphere. It follows from Eq. (1) that B z is 



time-independent; hence we do not show B z in the fig- 
ures. The geomagnetic field is set to Bq = 4.8 x 10~ 5 
Tesla, corresponding to an electron cyclotron frequency 
of 1.4 MHz, directed downward and tilted in the rrz-plane 
with an angle of 13 degrees (0.227 rad) to the z-axis, i.e., 
B = {B 0x ,B 0y ,B 0z ) = (sin0.227,0,-cos0.227)S . In 
our numerical simulation, we use 10 5 spatial grid points 
to resolve the plasma for < z < 400 km. The spa- 
tial derivatives are approximated with centered second- 
order difference approximations, and the time-stepping 
is performed with a leap-frog scheme with a time step of 
At = 8 x 10~ 9 s. 



III. NUMERICAL RESULTS 
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FIG. 2: The ion density profile, and the electric and magnetic 
field components at time t = 0.720 ms. The splitting of the 
wave is due to Faraday rotation. 

In the simulation, the EM pulse propagates without 
changing shape through the neutral atmosphere, until it 
reaches the ionospheric layer. At time t = 0.720 ms, 
shown in Fig. [2J the EM pulse has reached the lower 
part of the ionosphere. The initially linearly polarized 
EM wave undergoes Faraday rotation due to the differ- 
ent dispersion properties of the L-0 and R-X modes (we 
have adopted the notation "L-0 mode" and "R-X mode" 
for the two high-frequency EM modes, similarly as, e.g., 
Goertz and Strangeway [1995]) in the magnetized plasma, 
and the E y and B x components are excited. At t = 0.886 
ms, shown in Fig.[3j the L-0 and R-X mode pulses are in 
the vicinity of their respective turning points, the turning 
point of the L-0 mode being at a higher altitude than 
that of the R-X mode; see panel a) of Fig. [3J A closcup 
of this region, displayed in panel b), shows that the first 
maximum of the R-X mode is at z w 270.5 km, and the 
one of the L-0 mode is at z ~ 277 km. The maximum 
amplitude of the R-X mode is « 3 V/m while that of the 
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FIG. 3: a) The ion density profile, and the electric and mag- 
netic field components at time t = 0.886 ms. b) A closeup of 
the region of the turning points of the R-X and L-O modes. 
We see that the wave-energy of the L-O mode is concentrated 
into one single half-wave envelop at z ~ 277 km, while the 
turning point of the less localized R-X mode is at z ~ 270.5 
km. 



FIG. 4: a) The ion density profile, and the electric and mag- 
netic field components at time t = 0.948 ms. b) A closeup of 
the region of the turning points of the R-X and L-O modes. 
Here, the L-O mode oscillations at z fa 277 km are radiat- 
ing EM waves with perpendicular (to the z axis) electric field 
components. 



L-O mode is « 10 V/m; the latter amplitude maximum 
is in agreement with those obtained by Thide and Lund- 
borg, [1986], for a similar set of parameters as used here. 
The electric field components of the L-O mode, which at 
this stage are concentrated into a pulse with a single max- 
imum with a width of w 200 m, are primarily directed 
along the geomagnetic field lines, and hence only the E z 
and E x components are excited, while the magnetic field 
components of the L-O mode are very small. At t = 0.948 
ms, shown in panel a) of Fig. QJ both the R-X and L-O 
mode wave packets have widened in space, and the EM 
wave has started turning back towards lower altitudes. 



In the closeup of the EM wave in panel b) of Fig. dj one 
sees that the L-O mode oscillations at z w 277 km are 
now radiating EM waves with significant magnetic field 
components. Finally, shown in Fig. [5] at t = 1.752, the 
EM pulse has returned to the initial location at z = 50 
km. Due to the different reflection heights of the L-O and 
R-X modes, the leading (lower altitude) part of the pulse 
is primarily R-X mode polarized while its trailing (higher 
altitude) part is L-O mode polarized. In the center of the 
pulse, where we have a superposition of the R-X and L-O 
mode, the wave is almost linearly polarized with the elec- 
tric field along the y axis and the magnetic field along the 
x axis. The direction of the electric and magnetic fields 



here depends on the relative phase between the R-X and 
L-0 mode. 
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FIG. 5: The ion density profile, and the electric and magnetic 
field components at time t = 1.752 ms. 

In Fig. [6l panel a), we have plotted the electric field 
component E x at z — 270.50 km, near the turning point 
of the R-X mode and in panel b) we have plotted the E z 
component at z — 276.82 km, near the turning point of 
the L-0 mode. We see that the maximum amplitude of 
E x reaches 3 V/m at t = 0.87 ms, and that of E z reaches 
10 V/m at t — 0.9 ms. The electric field amplitude at 
2 = 270.50 km has two maxima, due to the L-0 mode 
part of the pulse, which is reflected at the higher alti- 
tude z = 276.82 km and passes twice over the altitude 
z = 270.50 km. We also observe weakly damped oscil- 
lations of E z at z = 276.82 km for times t > 1.05 ms, 
which decrease exponentially in time between t — 1.1 ms 
and t = 1.5 ms as E z oc exp(— jt) with 7 = 6.5 x 10 3 s _1 . 
We found from the numerical values that 7 w c/c„/2, 
where k n = dbxriio/dz ~ 4.6 x 10 _5 m _1 is the inverse 
ion density scale length at z — 277 km, but we are not 
certain how general this result is. No detectable mag- 
netic field fluctuations are associated with these weakly 
damped oscillations, and we interpret them as electro- 
static waves that have been produced by mode conver- 
sion of the L-0 mode. The amplitudes of the E x and E y 
components are also much weaker than that of the E z 
component for these oscillations. A closeup of these elec- 
trostatic oscillations at t = 1.152 ms is displayed in panel 

c) of Fig. [51 where we see that they have a wavelength 
of approximately 33 m (wavenumber 0.19 m _1 ). In panel 

d) of Fig. we have plotted the frequency / = uj/2tt as a 
function of the wavenumber k, where lu is obtained from 
the Appleton-Hartree dispersion relation [Stix, 1992] 




1 - 
0.5 - 
- 
-0.5 - 
-1 - 





d) 



276.5 


276.6 


276.7 


276.8 276.9 
z(km) 


277 


277.1 




R-X 










. 

L-0 












- ^ 










— e — 




Whistler 
/ 




\ 

Z 


X=33 m 




0.2 -0.15 


-0.1 


-0.05 


0.05 0.1 


0.15 


0.2 



k z (rrr 1 ) 

FIG. 6: a) The amplitude of the electric field component E x 
at z — 270.50 km, near the turning point of the R-X mode, 
and b) the amplitude of the electric field component E z at 
z = 276.82 km, near the turning point of the L-0 mode, c) A 
snapshot of low-amplitude electrostatic waves of wavelength 
A « 33 m (wavenumber k — 2-k/A ~ 0.19m -1 ), observed at 
time t = 1.152 ms, and d) Dispersion curves (lower panel) 
obtained from the Appleton-Hartree dispersion relation with 
parameters uj pc = 31.4 x 10 6 s _1 (5 MHz), uj cc = 8.80 x 10 6 s" 1 
(1.4 MHz) and 9 = 13° = 0.227 rad. We identify the high- 
frequency R-X and L-O modes, as well as the Z-mode which 
extends to the electrostatic Langmuir/Upper hybrid branch 
for large wavenumbers; the circles indicate the approximate 
locations on the dispersion curve for the electrostatic oscilla- 
tions shown in panel c). For completeness we also show the 
low- frequency electron whistler branch in panel d). 
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FIG. 7: The frequency spectrum (10-logarithmic scale) of a) 
the electric field component E x at the altitude z = 270.50 
km, and b) of E z at the altitude z = 276.82 km. 

Here A = [w c 2 e sin 4 6 + Alu- 2 (lu 2 - w 2 e ) 2 cos 2 9] 1 ' 2 , to pc 
(uj cc ) is the electron plasma (cyclotron) frequency, and 9 
is the angle between the geomagnetic field and the wave 
vector k, which in our case is directed along the z-axis, 
k = k z z. We use w pe = 31.4 x 10 6 s _1 (corresponding 
to / pc = 5 MHz), ui cc = 8.80 x 10 6 s^ 1 (correspond- 
ing to /ce = 1.4 MHz) and = 13° = 0.227 rad. The 
location of the electrostatic waves whose wavelength is 
approximately 33 m and frequency 5 MHz is indicated 
with circles in the diagram; they are on the same disper- 
sion surface as the Langmuir waves and the upper hybrid 
waves/slow Z mode waves with propagation parallel and 



perpendicular to the geomagnetic field lines, respectively. 
The mode conversion of the L-0 mode into electrostatic 
oscillations are relatively weak in our simulation of verti- 
cally incident EM waves, and theory shows that the most 
efficient linear mode conversion of the L-0 mode occurs 
at two angles of incidence in the magnetic meridian plane, 
given by, e.g., Eq. (17) in [Mj0lhus, 1990]. 

The nonlinear effects at the turning point of the L-0 
and R-X modes are investigated in Fig. [7] which displays 
the frequency spectrum of the electric field component E x 
at the altitude z — 270.5 km and of E z at the altitude 
z = 276.82 km. The spectrum shows the large-amplitude 
pump wave at 5 MHz and the relatively weak second har- 
monics of the pump wave at 10 MHz at both altitudes 
(the slight downshift is due to numerical errors produced 
by the difference approximations used in space and time) . 
Visible are also low-frequency oscillations (zeroth har- 
monic) due to the nonlinear down-shifting/mixing of the 
high-frequency wave field. 

IV. SUMMARY 

In conclusion, we have presented a full-scale numerical 
study of the propagation of an EM wave and its linear 
and nonlinear interactions with an ionospheric layer. We 
observe the reflection of the L-0 and R-X modes at dif- 
ferent altitudes, the mode conversion of the L-0 mode 
into electrostatic Langmuir/upper hybrid waves as well 
as nonlinear harmonic generation of the high-frequency 
waves. Second harmonic generation have been observed 
in ionospheric heating experiments [Derblom et al., 1989; 
Blagoveshchenskaya et al., 1998] and may be partially ex- 
plained by the cold plasma model presented here. 
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